function p=pickup(seismic_data,Tn,k)
%input:
%   seismic_data:地震数据
%   Tn:地震道集道数
%   k:地震道集序号
%output:
%   p:自动拾取地震初至时间矩阵
%email:2906310084@qq.com
%2020-6-4
    trace=seismic_data(:,(k-1)*Tn+1:k*Tn);%选取炮集记录数据
    for i=(k-1)*Tn+1:k*Tn
        p(i,1)=i;
        temp=1;
        while (temp>0&&temp<=2001)
            if(trace(temp,p(i,1)-(k-1)*Tn)~=0)%从1开始遍历，找到不为0的数，跳出循环
                break;
            end
            temp=temp+1;
        end
        t1=temp;%找到t1，起始
        while(temp>0)
            if(trace(temp,p(i,1)-(k-1)*Tn)*trace(temp+1,p(i,1)-(k-1)*Tn)<0)%从t1处接着遍历，找到零点位置
                break;
            end
            temp=temp+1;
        end
        t2=temp+1;%找到t2，终点
        plot(p(i,1),t1,'r*')
        plot(p(i,1),t2,'r*')%标记
        [~,idx]=min(trace(t1:t2,p(i,1)-(k-1)*Tn));
        p(i,2)=idx+t1;%得到波谷位置，并存入p矩阵
        plot(p(i,1),p(i,2),'rO')%标记波谷位置
    end
end
%     p(:,3)=floor(p(:,1)+0.5);
%     for i=1:size(p,1)%105
%         temp=floor(p(i,2));%取得时间
%              while(temp>0)
%                  if(trace(temp,p(i,3)-(k-1)*Tn)*trace(temp+1,p(i,3)-(k-1)*Tn)<=0)%判断是否找到位置
%                     break;
%                  end
%              temp=temp-1;
%              end
%         t1=temp-1;%找到时间 t1
%         temp=t1;
%             while(temp>0)
%                  if(trace(temp,p(i,3)-(k-1)*Tn)==0)%判断是否找到位置
%                     break;
%                  end
%              temp=temp-1;
%             end
%         t2=temp-1;
%         plot(p(i,3),t1,'r*')
%         plot(p(i,3),t2,'r*')
%         [temp,idx]=min(trace(t2:t1,p(i,3)-(k-1)*Tn));
%         t(i,1)=idx+t2;
%         plot(p(i,3),t(i,1),'rO')
%     end
% end
% end